Chaotic oscillation in attractive Bose-Einstein condensate under an impulsive force 
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For an attractive trapped Bose-Einstein condensate an imaginary three-body recombination loss 
term and an imaginary linear source term are usually included in the Gross-Pitaevskii (GP) equation 
for a proper account of dynamics. Under the action of an impulsive force, generated by suddenly 
changing the atomic interaction or the trapping potential, the solution of this complex GP equation 
for attractive interaction is found to lead to a very long term chaotic oscillation. 
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Since the successful detection [jy, |2j of Bose-Einstein 
condensates (BEC) in dilute bosonic atoms employing 
magnetic trap at ultra-low temperature, there has been 
great theoretical and experimental interest in conden- 
sates with attractive interaction [Q. For attractive in- 
teraction the condensate is stable for a maximum critical 
number Acr of atoms When the number of atoms 
increases beyond this critical number, due to interatomic 
attraction the condensate collapses emitting atoms until 
the number of atoms is reduced below Acr and a stable 
configuration is reached. With a supply of atoms from 
an external source the condensate can grow again and 
thus a series of collapses can take place, which was ob- 
served experimentally in the BEC of 7 Li with attractive 
interaction Theoretical analyses based on the mean- 
field Gross-Pitaevskii (GP) equation || Q also confirm 
the collapse |,|,§. 

In a recent classic experiment by the JILA group on 
the BEC of 85 Rb, a sustained harmonic oscillation of an 
attractive condensate had been observed and measured 
when the scattering length of atomic interaction was sud- 
denly changed (?]]. It is possible to manipulate the scat- 
tering length by an external magnetic field via a Feshbach 
resonance ||, |[ . The sudden change in scattering length 
or the harmonic oscillator trapping potential constitute a 
general class of impulsive force in the GP equation which 
can be realized experimentally. We show that under the 
action of an impulsive force the solution of the GP equa- 
tion for attractive interaction with a pentic loss term and 
a linear source term can exhibit sustained chaotic oscil- 
lation. The imaginary pentic term in the GP equation 
accounts for three-body recombination loss and the imag- 
inary linear term accounts for growth due to an external 
source of atoms. 

There have been several studies of chaotic dynamics 
in the nonlinear Schrodinger p0[ , Landau-Ginzburg p"T| , 
and the GP |TJ, [TJ) equation. But all these studies 
were motivated from a numerical or mathematical point 
of view without phenomenological consequence. The 
chaotic oscillation we study here may be observed in 
the laboratory by slightly manipulating the strength of 
the linear source term responsible for supplying atoms 
to the condensate. In the absence of any source term 
or for a weak source term the oscillation is periodic af- 
ter the application of the impulsive force. With a small 



increase in the source term the motion turns out to be 
chaotic. The JILA group has performed a classic experi- 
ment where they have been able to observe and measure 
the frequency of oscillation of a condensate after the ap- 
plication of an impulsive force generated by changing the 
scattering length to a negative value ]tJ . Such an experi- 
ment might be able to detect the transition from periodic 
to chaotic oscillation of a condensate under the action of 
an impulsive force. 

We consider for our present study on the numerical 
solution [|l4|, [l5| of the GP equation || for a spherically 
symmetric harmonic trap. In the GP equation we in- 
clude a pentic three-body nonlinear recombination term 
which accounts for the decay of the strongly attractive 
condensate and a linear source term. 

The GP equation in this case in dimensionlcss units 
can be written as |^|, || 
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ip(x,t) 
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tp(x,t) = 0, 



(1) 



where n = Na/l, 7 is the coefficient of the linear source 
term and £ is the coefficient of the pentic three-body re- 
combination term. Here the distance x, time t, and the 
spherically symmetric wave function (f>(x,t) = (p(x,t)/x 
are expressed in units of l/s/2, u^ 1 , and (VStt^ 3 ) -1 , 
where I = (h/muj) and u> is the frequency of the har- 
monic oscillator trap, m is the mass of a single atom, A^ is 
the number of atoms in the condensate, and a the scat- 
tering length of atomic interaction. The normalization 
condition for the wave function is 



dx\cp(x,t)\ 2 = 1. 



(2) 



We solve Eq. (Q) above by the Crank-Nicholson time 
propagation Q after discretization with 7 7^ and £ 7^ 
starting from the known harmonic oscillator solution for 
n = £ = 7 = 0. The x discretization was performed 
with space step 0.1 up to a maximum x of 40. The time 
step was taken to be 0.01. The nonlinearity constant n 
was increased by steps of 0.0001 until the desired value 
is reached. Once the final nonlinearity is reached the so- 
lution is then stabilized by iterating 100000 times which 
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corresponds to an interval of time t — 1000. This removes 
any transient behavior in the solution which is important 
for a study of chaos. Then the wave function is prepared 
for the simulation of chaos after the application of the 
impulsive force. The oscillation of the system after the 
application of the impulsive force is best studied by con- 
sidering the time evolution of the root mean square (rms) 
radius X and its time derivative X. 

For attractive interaction a stable condensate can be 
formed for the minimum value —0.575 for n ji], |[ rc|. We 
consider the following three numerical simulations for the 
study of chaos: (i) On a preformed attractive condensate 
with n — —0.4 we suddenly double the harmonic oscil- 
lator term from x 2 /4 to x 2 /2 and study the resultant 
oscillation of the system for different 7. On a preformed 
repulsive condensate with n = +0.4 we suddenly change 
the sign of the scattering length a to (ii) —a and also to 
(iii) —2a and study the resultant oscillation of the system 
for different 7. In all simulations we take £ = 0.0004. In 
cases (i) and (ii) the final value of n(— —0.4 > —0.575) 
permits stable condensate and in case (iii) final value of 
ri(= —0.8 < —0.575) does not permit a stable solution of 
Eq. (|l|). However, the experiment of JILA |?J has shown 
that in case (iii) the number of particles can not decay 
immediately below the stability limit and there could be 
sustained oscillation of the system with a value of n less 
than —0.575. We find that depending on the value of 
the source term 7 there could be pronounced chaos in all 
three cases above. 

Once the wave function is prepared for the study of 
chaos as described above we inflict the change corre- 
sponding to case (i), first, with 7 = 0. Only periodic os- 
cillation of rms radius X and its derivative X is observed. 
With a slight increase in 7 to 0.0002, chaotic oscillation 
is obtained. The change from periodic to chaotic motion 
is best illustrated by plotting X versus X. For a periodic 
motion a closed loop appears in the phase space plot of 
X versus X, whereas for a chaotic motion a strange at- 
tractor appears in the phase space plot. This is shown in 
Figs. 1 (a) and (b) for 7 = 0.0001 (periodic) and 0.0002 
(chaotic) , respectively, for an interval of simulation time 
t = 200 after inflicting the impulsive force. With the 
increase of 7, the periodic motion of Fig. 1(a) changes 
to chaotic motion in Fig. 1(b). The present time is ex- 
pressed in units of u^ 1 . In a typical experimental situ- 
ation the harmonic oscillator trapping frequency v ~ 50 
s , and hence lo^ 1 = (2ttu)~ 1 ~ 0.005 s. Thus the sim- 
ulation time of t = 200 corresponds to 1 s which is inside 
the experimental observation period of a typical set up 
0. We continued the simulation till t — 20000 and the 
robust chaotic attractor seems to stay for ever, although 
it moves slowly to a smaller value of X with time. This 
is exhibited in Fig. 1 (c) where we plot X versus X for 
< t < 20000 which corresponds to an interval of 100 s. 

Next we perform a similar simulation as in Fig. 1 
above for the case (ii) where the sign of the scattering 
length is suddenly changed from positive to negative ex- 




FIG. 1: X versus X in case (i) when the harmonic oscillator 
potential is suddenly doubled in Eq. (|l|) for £ = 0.0004, n = 
-0.4 and (a) 7 = 0.0001, t = 200, (b) 7 = 0.0002, t = 200, 
and (c) 7 = 0.0002, t = 20000. 



plotting a Fcshbach resonance as in the experiment of 
JILA H. In simulation this corresponds to changing 
n from 0.4 to —0.4 (repulsive to attractive). The fi- 
nal n(= —0.4 > —0.575) allows for a stable attractive 
condensate to be formed. Again for 7 = and 0.0001 
periodic motions are obtained which change to chaotic 
motion as 7 is increased. In Fig. 2 we plot X versus X 
for 7 = 0.0001 and 7 = 0.0002 as in Fig. 1. In Fig. 2 (a) 
we show the periodic motion for 7 = 0.0001 for time t up 
to 200. Fully chaotic motion is obtained for 7 = 0.0002. 
The chaotic motion for the initial interval of time 200 is 
shown in Fig. 2 (b) and that for an interval of 20000 
units of time is shown in Fig. 2 (c). The robust chaotic 
attractor has started to develop in Fig. 2 (b) and its 
long-term evolution is shown in Fig. 2 (c). 

Finally, we consider the case (iii) above. In this case 
the nonlinear term n has been suddenly changed from 
0.4 to —0.8 (repulsive to attractive). However, in this 
case the final n(= —0.8 < —0.575) does not allow a sta- 
ble condensate to be formed. Such a strongly attractive 
condensate has been created and observed in the labo- 
ratory in the experiment of JILA J7J. In this case the 
condensate is unstable and due to interatomic attraction 
it starts to shrink in size or collapse after the sudden 
change in the nonlinear term. Consequently, as the cen- 
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FIG. 2: X versus X in case (ii) when the sign of the scattering 
length a was suddenly made negative from positive in Eq. (fj ) 
for £ = 0.0004, n = 0.4 and (a) 7 = 0.0001, t = 200, (b) 
7 = 0.0002, t = 200, and (c) 7 = 0.0002, t = 20000. 



FIG. 3: X versus X in case (iii) when the scattering length 
a was suddenly changed to —2a in Eq. ([j]) for £ = 0.0004, 
n = 0.4 and (a) 7 = 0, t = 200, (b) 7 = 0.0002, i = 200, and 
(c) 7 = 0.0002, t = 20000. 



tral density of the condensate increases, it starts to emit 
particle through small explosions and tries to attain a 
more stable configuration with a smaller number of par- 
ticles [jjj . The condensate exhibits oscillation during this 
process of collapse and explosion but it may need a large 
amount of time before attaining the critical size with A^r 
atoms. No periodic oscillation was observed in this case 
even for 7 = 0. In Fig. 3 (a) we plot X versus X for the 
first 200 units of time for 7 = 0, which shows a chaotic at- 
tractor, which continues to exist for any positive nonzero 
7. In Fig. 3 (b) we plot the same for 7 = 0.0002, which 
shows the strange attractor for the first 200 units of time. 
In Fig. 3 (c) we exhibit the long-term behavior of this 
attractor for 20000 units of time. In the long term the 
attractor is more spread in this case compared to Figs. 
1 (c) and 2 (c). In case (iii) in the long term the chaotic 
attractor moves to large X whereas in cases (i) and (ii) 
it moves to smaller X. 

The presence of chaos is characterized by finding the 
Lyapunov exponents [ fl6| . Since the quantities X and X 
are not direct solutions of the nonlinear GP equation (|l|) 
and are essentially the expectation values, we use time 
series analysis to calculate the Lyapunov exponents jl6) . 
In all the cases, we use the velocity-variable (X) data 
collected with the time interval St = 0.1. The reason to 



choose X is that it oscillates with time only around the 
steady average zero value, whereas the time average of 
X is not stationary and changes with time. From one set 
of data for X(t) we construct other independent set(s) 
by allowing a time lag r. The number of such indepen- 
dent sets constitute the embedding dimension m of the 
data sets for X(t) to be analyzed. The optimal time lag 
is found from the estimate of mutual information which 
shows the independent nature of the reconstructed data 
sets . After some experimentation we find that m = 3 
leads to reasonable values for all the Lyapunov expo- 
nents. All the calculations of the Lyapunov exponents 
are performed with m = 3 using the algorithm by Sano 
and Sawada Q. 

The calculation leads to three exponents in each of the 
cases (i), (ii), and (iii) with 7 = 0.0002 shown in Figs. 1 
(c), 2 (c), and 3 (c), respectively. In the case of (i), the 
time lag r in the calculation of Lyapunov exponents was 
found to be 0.6 and in cases (ii), and (iii) it is found to 
be 0.9. In Fig. 4 we plot the largest Lyapunov exponent 
versus the time used for the calculation in the three cases 
above. In all three cases the largest exponent is found to 
lead to a convergent positive value at large finite time 
which shows the existence of chaos in all the cases . We 
also found that the oscillation in Fig. 3 (a) corresponding 
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to case (iii) with 7 = is also chaotic. 
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FIG. 4: The Lyapunov exponents for the chaotic attractors 
of Figs. 1 (c), 2 (c), and 3 (c) labeled (i), (ii), and (iii), 
respectively. 

There has been confirmation of chaos in the nonlinear 
Schrodinger equation ]ic| and in the Landau-Ginzburg 
equation |pd| with similar nonlinear terms as in the GP 
equation, although the details of the two equations are 
different. Also, there has been theoretical prediction of 
chaos in a coupled set of GP equations Jl^] . The numer- 
ical study of chaos in a collapsing Bose-condensed gas 
in Ref. fig ] is worth mentioning. In that work chaos 
was confirmed in the normal time evolution of the GP 
equation (|l|) without any external impulsive force. The 
chaotic attractor in that study was found to stay only 
for a short interval of time (~ 500) , where in the present 



study with impulsive forces the chaotic oscillations are 
found to appear immediately after the application of the 
impulsive force and stay for more than 20000 units of 
time t. This makes the analysis of the chaos via Lya- 
punov exponent more reliable. 

The appearance of chaos in nonlinear dynamics is 
of interest from a theoretical point of view. Here we 
have demonstrated the chaotic oscillation in an attrac- 
tive BEC under an impulsive force using the mean-field 
GP equation. The use of the GP equation in this study 
is justified as this equation produces a faithful represen- 
tation of the BEC for both repulsive and attractive in- 
teractions ||, |(| . We also performed simulation when the 
condensate is repulsive after the application of the im- 
pulsive force. No chaotic oscillation was detected in that 
case. From this it seems that the final-state interatomic 
attraction plays an important role in the generation of 
this chaotic dynamics. 

In summary, from a numerical simulation based on the 
solution of the Gross-Pitaevskii equation (jl]) for attrac- 
tive interaction with an absorptive pentic three-body re- 
combination term (£) and a linear source term 7, we find 
that sustained chaotic oscillation can result in a BEC 
under the action of an impulsive force generated by sud- 
denly changing the interatomic scattering length or the 
harmonic oscillator trapping potential. 
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